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Abstract 

With a nonequilibrium relaxation method, we calculate the dynamic 
critical exponent z of the two-dimensional Ising model for the Swendsen- 
Wang and Wolff algorithms. We examine dynamic relaxation processes 
following a quench from a disordered or an ordered initial state to the 
critical temperature T c , and measure the exponential relaxation time of 
the system energy. For the Swendsen-Wang algorithm with an ordered 
or a disordered initial state, and for the Wolff algorithm with an ordered 
initial state, the exponential relaxation time fits well to a logarithmic 
size dependence up to a lattice size L = 8192. For the Wolff algorithm 
with a disordered initial state, we obtain an effective dynamic exponent 
Zcxp = 1.19(2) up to L = 2048. For comparison, we also compute the 
effective dynamic exponents through the integrated correlation times. In 
addition, an exact result of the Swendsen-Wang dynamic spectrum of a 
one-dimension Ising chain is derived. 

1 Introduction 

In recent two decades, cluster algorithms have played an important role in sta- 
tistical physics due to their reduced critical slowing down, improved computa- 
tional efficiency, and interesting dynamical properties. Among these dynamical 
properties, the dynamic critical exponent z is the center of attraction, which 
describes the divergent correlation time. 

There are various ways to calculate the dynamic critical exponent z, for 
example, through the exponential decay of the time correlation of a finite system 
in equilibrium or from the dynamic scaling behavior in nonequilibrium 
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states [3JElini ■ In calculating the time correlation in equilibrium, the difficulty is 
that one can hardly reach a very large lattice. The advantage for computing the 
dynamic exponent from a nonequilibrium relaxation process is that the finite 
size effect is more or less negligible, since the spatial correlation length is small 
in the early stage of the dynamic relaxation. Such a nonequilibrium approach, 
however, becomes subtle for the cluster algorithms, for the dynamic exponent 
z is believed to be close to zero. In addition, it is also somewhat controversial 
in defining a Monte Carlo time for the Wolff algorithm. 

In a recent article & , an attempt is made to estimate the dynamic exponent z 
of the Wolff algorithm from the finite size scaling behavior in a nonequilibrium 
state. A vanishing z value is claimed. To our understanding, however, the 
identification of the dynamic scaling behavior there seems to be not appropriate. 
On the other hand, although the cluster algorithms are known to be very efficient 
in reducing critical slowing down with a small dynamic exponent z, it has not 
been rigorously studied what precise values the dynamic exponent z takes for 
different variants of the algorithms. This is important in theory and application 
of the cluster algorithms. 

In this paper, we will calculate the dynamic exponent z for the Wolff an d 
Swendsen-Wang [H| algorithms, respectively, through the exponential relaxation 
time and integrated correlation time [§] of the system energy in nonequilibrium 
relaxation processes. Compared with methods based on computations of time 
correlation functions in equilibrium, much larger system sizes can be reached 
in the nonequilibrium dynamic approach, especially in the case of the two- 
dimensional Ising model, where the system energy in the equilibrium state is 
known exactly. Compared with the methods in Ref. jSJ, the system energy in 
our calculations is self-averaged and thus much less fluctuating in simulations 
of large lattices. 

The paper is organized as follows. In Sec. 2, the general theory of the 
spectrum of the Monte Carlo dynamics is described, and in Sec. 3, an exact 
calculation of the spectrum of the one-dimensional Ising chain is formulated. In 
Sec. 4, simulation setup is discussed in detail, and in Sec. 5, numerical results 
are presented. 



2 Spectrum of Monte Carlo Dynamics 

In order to justify our method, we first look at the spectrum of a Markov chain 
Monte Carlo dynamics JUj and its relation to observables, i.e., the equilibrium 
and nonequilibrium relaxation functions. Let W be a transition matrix of an 
irreducible, aperiodic, and reversible Markov chain with an equilibrium (invari- 
ant) probability distribution p. We have a detailed balance equation between 
W and p, 

p i W ij =p j W jt . (1) 
This equation implies that the following matrix is symmetric: 

Si^p^WvpJ 1 ' 2 , (2) 
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thus, the eigenvalues A m of S are real. Due to the conservation of the total 
probability, it can also be shown that |A m | < 1. Let the eigenvectors of S be 
Ui m for eigenvalue A m , then the left and right eigenvectors of W is Xi — p\^ 2 Ui m 
and y-i = p\ 1 ^ 2 Ui m , respectively, such that 

xW = X m x, Wy = X m y. (3) 

The equilibrium distribution corresponds to Ao = 1, x^ = p, and yf 3 ^ = 1. The 
next eigenvalue Ai nearest to 1 controls the rate of convergence. We define the 
exponential relaxation time r (in units of one Monte Carlo step or attempt) by 
Ai = exp(-l/r). 

We can represent the relaxation of a general observable Q in terms of the 
initial distribution p(0) or equilibrium distribution p and the eigenspectrum of 
S as 

(Q(*)> P (0) = ^XldkCk, (4) 
k 

{Q(t)Q(0)) eq = J2^l (5) 
k 

where the averages are over the initial distribution and equilibrium distribution, 
respectively, and 

Ck = ^2Pi /2 QiUik, (6) 

i 

dk = ^2p7 1/2 Pt(®) u ik- (7) 

i 

We define the normalized relaxation function f(t) to be linear in (Q(t)) or 
(Q(i)Q(O)} such that /(0) = 1 and /(oo) = 0, e.g., /(t) = «Q(t)>-(Q(oo)»/«Q(0)>- 
(<3(oo))). The integrated correlation time is defined as 



Tint 



£/(*)• (8) 



We note that the integrated correlation time depends not only on the observ- 
able Q but also on the dynamics and the full eigen spectrum. The integrated 
correlation time for the equilibrium correlation and noncquilibrium relaxation 
is not the same. On the other hand, the exponential relaxation time, defined in 
the large time limit, f(t) ~ exp(— t/r), is an intrinsic property of the Markov 
chain, ft is the same for both the equilibrium and nonequilibrium situations. 

3 Exact Calculation in One Dimension 

It is instructive to look at the eigen spectra of the Swendsen-Wang and Wolff 
dynamics in one dimension (ID). We consider the Swendsen-Wang dynamics in 
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ID with an open boundary condition. Consider spin <ii at a one-dimensional 
lattice site i = 1, 2, • • •, L, L + 1. We define the link variable bi = 1 — 8 ai _ ai+1 . 
The energy of the system is 

L L 

£(cr) = -J^o-jCTj+i = 2 J^b* + const. (9) 

i=l i=l 

We introduce the bond variables rii = 0, 1 representing the absence or presence of 
a bond in the Swendsen-Wang dynamics, then the joint probability distribution 
of both the spin and bond is proportional to 

L 

P{a,n) = n^^+i^.i + (1 -p)S ni ,o], (10) 
i=i 

where p = l-exp[-2J/(fc B T)]. The marginal distribution of the spins is given 
by J^n p ( a i n ) = Pi' 7 ) = EL exp[-2 Jbi/(k B T)]. The distribution of the bonds 
is the special case of the Fortuin-Kasteleyn formula, P{n) = P(a, n) — 
p N "{\ - p )L-N b2 L-N b +i^ Nb = £\ ni . The Swendsen-Wang algorithm is to 
apply alternatively two conditional probabilities, sample bonds given the spins 
P(n\a) = P(a, n)/P(a), and sample spins given the bonds P(a\n) = P(a, n)/P(n). 
Instead of using the spin variables, it is more convenient to use the link variable 
bi. In terms of bi, the conditional probabilities are simple products: 

P(n|6) = n/(6 i> n<) J /= ( 1 \ P P Q ) , (11) 

P(&|n) = n<?KM: 9=( 1 { 2 (12) 

The transition probability from a given link configuration to another link con- 
figuration is 

W(b^b')=Y,P(b'\n)P(n\b)=l[ Wbt ^ ™=(% (13) 

n i \ / / / 

The matrix w has eigenvalues 1 and p/2, with left eigenvectors v^> =(1,1 — p) 
and = (1,-1), respectively. We note that the full matrix W for the whole 
system is a direct product from the contribution of each site. Thus, the eigen 
values of W are A m = (p/2) m |llj . with L\/(m\(L — to)!) fold degenerate eigen 
vectors, \\i v^ ki '(bi), ki — 1 or 2, with L — m terms of choices for fcj = 1 and to 
terms for ki — 2, where m = 0, 1, 2, • • • , L. The eigenvalue Ao = 1 corresponds 
to the equilibrium state with a left eigenvector P{u). The next eigenvalue 
Xi = p/2 = exp(— 1/r) gives the relaxation time. We note that rest of the decay 
times t/to are well-spaced. 

It is easy to write down the probability distribution of the domain length 
since each link evolves independently. Let Pi (t) be the probability for observing 
a domain of a length I with — spins terminated by + spins, then 

P l {t)=p(b=\,t)[\-p{b=l,t)] l -\ (14) 
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Figure 1 : Eigenvalues of the transition matrix of the Swendsen-Wang dynamics 
in the ID Ising chain with the open boundary condition (solid lines) and the 
periodic boundary condition (dotted line with circles) for a lattice size L = 8. 

where p(t) = p(0)u/ is the probability of the link variable takes the value 1 at 
step t. Similar result using a continuous time dynamics is given in Ref. |12| . 

When a periodic boundary condition is used, we are no longer able to find 
the eigen spectrum analytically. In Fig. ^ we show the numerical results by di- 
agonalizing the transition matrix on an L = 8 chain with the periodic boundary 
condition and compare with the open boundary condition result. We can make 
several interesting observations. The eigenvalue Ai = p/2 is always present for 
both periodic and open boundary conditions for any lattice size L. However, 
we are unable to prove it rigorously. With the periodic boundary condition, the 
eigen value p/2 no longer corresponds to the slowest mode. It seems reason- 
able that the simple spectrum (p/2) m is a good approximation if the correlation 
length £ ~ exp(2 J/(fcsT)) is much smaller than the system size L, thus it is 
the correct spectrum in the thermodynamic limit. For finite sizes when £ is 
comparable to size L, the degenerate spectrum splits and rejoins at T = Q. 

4 Simulation Setup 

The two-dimensional (2D) Ising model is prepared initially in a random state 
with a zero magnetization or at the ground state, but evolves at the critical 
temperature. After a certain time before reaching equilibrium, one can observe 
an exponential decay of the system energy, which satisfies the following equation: 



E{t) « Ae~ t/T +E{oo), 



(15) 
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where r is the so-called exponential relaxation time. The exponential relaxation 
time r is an intrinsic property of the Monte Carlo algorithm, which is dehned by 
the first excited eigenvalue Ai of the transition matrix, and should be indepen- 
dent of the initial states in the simulations. At the transition temperature, the 
dynamic scaling theory predicts that r diverges according to r ~ L Zcxp in the 
thermodynamic limit. This defines the exponential dynamic critical exponent 
z cxp . On the other hand, the integrated correlation time is defined as 



t= 



Similarly, one may define the integrated dynamic critical exponent z- m t through 
Tint ~ L Zint . In order to get a more accurate result, we use the exact value of 
E(oo) for the 2D Ising model |T5] . 

The Wolff algorithm exhibits an important difference compared to other al- 
gorithms. It only updates the spins belonging to a certain cluster around the 
seed spin at each Monte Carlo step, while other algorithms sweep the whole 
lattice. For a fair comparison with the Swendsen-Wang or single-spin-flip algo- 
rithms, we need to rescale the Wolff Monte Carlo steps. Specifically, the Monte 
Carlo time t' in the Wolff algorithm should be transformed to t 

-fsp (it) 

t"=i 

where t' or t" is the Monte Carlo time step of the Wolff single cluster flip (i.e. 
the number of clusters flipped so far), C(t") is the average size of the cluster 
at step t", and L d (d = 2) is the total number of spins of the system, t is 
proportional to the actual CPU time. This newly scaled time t should then be 
used in Eq. (|15|) for the exponential relaxation time of the Wolff dynamics. 
The integrated correlation time should also be changed to 

T ^ - E (0) - E(oo) X L d ■ (18j 

In the Swendsen-Wang algorithm, there is no complication in the definition 
of time, we straightforwardly use Eq. 1)15(1 and (|16JI to calculate the dynamic 
exponent z cxp and z- mt . 

Our definition of the transformed time t is slightly different from that in 
Ref. 0. In Ref. 0, the time t is defined by t — t'C(t'). If C(t') takes a scaling 
form C(t') cx t' a , two definitions coincide. However, the disadvantage of the 
definition in |H] is that one needs to assume or know the scaling behavior of the 
cluster size C(t ). Our definition of the time t has clear physical meaning, i.e., 
whenever the number of the flipped spins reaches L d , it counts as a unit time. 
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5 Results 



We use the standard Hamiltonian of the two-dimensional Ising model, 

- /3H = K^2,a i a j . (19) 

m 

Here (3 = l/(ksT) and K = J/(fc,gT), ks is the Boltzmann constant, T is 
temperature and J is the interaction energy between two spins. Spins <Ji take 
only values +1 and —1. The site i is on a square lattice with periodic boundary 
conditions. 

Let us start our numerical simulations with the Swendsen-Wang algorithm. 
To calculate the dynamic exponent z cxp of the Swendsen-Wang algorithm, we 
have used lattices L = 64, 128, 256, 512, 1024, 2048, 4096, 8192, and each 
has 2 24 , 2 21 , 2 19 , 2 17 , 2 15 , 5 x 10 4 , 4 x 10 4 , 2 x 10 4 runs respectively, and the 
maximum Monte Carlo time steps of each lattice are 60, 70, 70, 80, 90, 100, 
100, 110. Here two different initial temperatures T — oo and T — are used 
in the simulations. The system evolves at the critical temperature. From the 
exponential decay of the system energy in Eq. (| 1 51) , one measures the relaxation 
time r. The results are plotted on a linear-log scale in Fig. [5J The relaxation 
times are almost identical for both the hot start and cold start. Obviously, an 
approximate linear behavior is observed for the lattice size L > 100 in Fig. [21 
In other words, a logarithmic dependence r ~ \nL gives a better fit to the 
numerical data. If we take into account the data of smaller lattices, the In- 
dependence of the relaxation time could be given by r ~ (lnL) 12 . 

If we analyze the data in the form of a power-law dependence, we found that 
the effective exponent z oxp decreases continuously with increasing lattice sizes, 
reaching 0.18 at the largest size simulated. This is shown with a log-log scale in 
Fig. |3| and it strongly suggests that the relaxation time does not follow a power 
law with a small but finite exponent. Our extensive data with large system 
sizes thus agree with the conclusion of Heermann and Burkitt |14|. Previous 
calculations [SJ ED EH E3 EH EH g ave various values ranging from 0.35 to 0.2. 
This appears to be an effect of finite lattice sizes. 

For comparison, we also calculated the integrated relaxation time for the 
Swendsen-Wang dynamics, it is nearly a constant around r; n t « 3.1. This 
implies that z; nt ss 0, it is consistent with z cxp . 

We now investigate the Wolff dynamics with the fully ordered state at T = 
as the initial state, and evolve the system at the critical temperature. The 
lattice sizes are L = 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192 with 
2 2 9j 2 27 ; 2 m 2 2i ; 2 2 0j 2 20 ; 2 u 5 x 1q5j 2 x 1Q 5 ) g x 2 x 10 4 independent 

runs, respectively. The maximum Monte Carlo time steps of each lattice are 80, 
100, 150,180, 200, 200, 200, 220, 250, 250, 250. We observe that the dynamic 
behavior here is very similar to that of the Swendsen-Wang dynamics. As shown 
in Fig. 01 the correlation time exhibits a logarithmic size dependence even from a 
relatively small lattice size. The integrated relaxation time is nearly a constant, 
fiat ~ 1-17 when the lattice size is larger than 256. This implies that Zj nt « 0. 
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Figure 2: r eX p versus L in linear-log scale for the 2D Ising model with the 
Swendsen-Wang dynamics. The diamonds are for random initial configurations, 
and circles are for the ordered initial state. A straight-line fit gives r oxp = 
-3.42 +1.98 In L. 
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Figure 3: T cxp versus L in double logarithmic scale for the 2D Ising model with 
the Swendsen-Wang dynamics. 
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Figure 4: r versus L in linear-log scale for the 2D Ising model with the Wolff 
dynamics starting from an ordered initial state. The diamonds are exponential 
relaxation times, and a straight-line fit gives T cxp = 1-09 + 0.85 InL. The circles 
are integrated relaxation times. They reach a steady value near 1.17 for lattice 
sizes larger than 256. 

It is intuitively understandable that the Wolff dynamics starting from a fully 
ordered state is similar to the Swendsen-Wang dynamics. When the Swendsen- 
Wang algorithm is initialized in the order state (T = 0), and then evolves at 
critical temperature, it forms very large clusters which dominate the evolution 
of every Monte Carlo step, while other small clusters' effect can be neglected. 
Thus its behavior is much like a Wolff algorithm initializing at the ordered state, 
for it also has a very large cluster dominating the dynamic evolution. For the 
Swendsen-Wang algorithm with a disordered initial state, our simulations show 
that the clusters grow rapidly, and the relaxation time is almost the same as 
that with an ordered initial state. 

For Wolff dynamics with completely disordered initial states, we have used 
lattices L = 64, 128, 256, 512, 1024, 2048 with 2 24 , 2 23 , 2 20 , 2 17 , 2 15 , 2 8 
runs respectively. The maximum Monte Carlo time steps of each lattice are 
450, 1337, 4770, 18000, 69910, 300000 in the original time unit of t'. For the 
Wolff algorithm with a disordered initial state, the dynamic behavior is rather 
complicated. In order to have a better comparison of the ordered and disordered 
initial starts, we present the dynamic relaxation of the system energy for the 
lattice size L = 2048 with two different initial conditions in Fig. 03 The figure 
shows that the dynamic relaxation with a disordered start is much slower. 

In Fig. [13 T cxp and Tj nt are plotted as functions of the lattice size L in log-log 
scale. For r exp , an approximate power-law behavior is observed, and from large 
lattice sizes one derives a dynamic exponent z cxp rs 1.0. To obtain a better 




Figure 5: \E(t) — E(oo)\ versus t in log-linear scale for the 2D Ising model with 
transformed t for the Wolff dynamics. The lattice size is L = 2048. 



value of z, one may consider corrections to scaling. For example, assuming 
T CX p ~ L z (l + c/L s ), the fitted dynamic exponent is z cxp = 1.19(2). Similarly we 
estimate z- m t = 0.29(3) from Tint- z cxp is much larger than the results obtained 
in equilibrium correlation functions |15j . Both of these values are quite different 
from that reported in Ref. [Sj , where it is concluded that the exponent z is very 
close to zero. 

To understand the difference between our analysis and that in Ref. [Hj, we 
have also performed the scaling plot of Fig. 1 in Ref. 6 with our numerical 
data. Indeed, the scaling collapse is observed for large lattices. According to 
our scaling analysis, however, a dynamic exponent z ~ 1.7 should be extracted. 
Following the definition in Ref. t = t' (C(t'))/L d with t' being the Monte 
Carlo time step of the Wolff single cluster flip, and then r = T'(C)/L d . As- 
suming that (C) behaves like susceptibility, i.e., (C) ~ L 1 ^, one may deduce 
z = z '-(d-~f/v) = z'+2(Y H -d). In Ref. 0, this is written as z = z'-(2Y H -d). 
Together with other inconsistent formulations, a dynamic exponent z ~ is de- 
rived. 

In numerical simulations in equilibrium, one tends to conclude that the dy- 
namic exponent z of the Wolff algorithm is close to zero. This is in agreement 
with our results from the dynamic relaxation starting from an ordered initial 
state. Compared with the dynamic relaxation of the Wolff algorithm with an 
ordered initial state, why does the dynamic relaxation with a disordered initial 
state show an anomalous behavior? Our conjecture is that the eigenvalues A m of 
the transition matrix of the Wolff algorithm are rather dense. Within a rather 
long time, the contribution of the higher eigenvalues to the dynamic observable 
will not be suppressed in the dynamic relaxation starting from a disordered 
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Figure 6: r versus L in double logarithmic scale for the 2D Ising model with the 
Wolff algorithm starting from disordered states. The diamonds are exponential 
relaxation times r cxp . The circles are integrated relaxation times. 

state. The dynamic exponent measured in this paper and in Ref. [H] is actually 
an effective one. This also explains why the dynamic exponent extracted from 
the exponential decay of the system energy is somewhat smaller than that from 
the dynamic scaling behavior in the relatively short time regime in Ref. [H]. 
Our numerical simulations and data analysis actually show that if possible, one 
should avoid starting the simulations with the Wolff algorithm from a disordered 
state. 

Finally, we should mention that for the dynamic relaxation of the Wolff 
algorithm with a disordered initial state, it is hard to measure the relaxation 
time T exp in an extremely long time regime where the contribution of the higher 
eigenvalues is already suppressed, since \E(t) — E(oo)\ is too much fluctuating. 

6 Conclusion 

In summary, we have computed both the exponential and integrated relaxation 
times for the Wolff and Swendsen-Wang algorithms, with both disordered and 
ordered initial states. For the Swendsen-Wang dynamics, the exponential re- 
laxation time shows a logarithmic dependence on the lattice size L for both 
initial states, while the integrated relaxation time tends to a constant. This 
is a strong evidence that z cxp = for the Swendsen-Wang algorithm. For the 
Wolff dynamics with an ordered initial state, the results are similar to those 
of the Swendsen-Wang dynamics. For the Wolff dynamics with a disordered 
initial state, however, the dynamic relaxation is very slow, and it takes a very 
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long time to approach equilibrium. If one measures the relaxation time r exp 
in a reasonable time regime in the simulations, an effective dynamic exponent 
z C xp — 1-19(2) is obtained. 
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